### Load libraries
library("foreign")
library(car)
library(MASS)

### Clear space
rm(list = ls())

### Load data from the internet
#### First download corrected data
setwd("")
df <- read.dta("bjps_powell_tucker_dataset.dta")

### Limit data to their sample of 89 elections
df.fh <- df[df$EEIndicator=="1",]

#######################
#### Make Figure 1 ####
#######################

### Generate labels
df.fh$pclabel <- paste(df.fh$Country, df.fh$Year2, sep=", ")

### Bivariate plot
par(mar=c(5, 5, 1, 2), oma=c(0,0,0,0), cex.main=0.9, cex.axis=0.9, cex.lab=1)
plot(df.fh$GDPChange1989, df.fh$TypeAVol, xlab="GDP Change From 1989", ylab="Replacement Volatility",main="", pch=ifelse(df.fh$Country=="Bosnia-Herzegovina", 19, 1), bty="n", ylim=c(0,90), col=ifelse(df.fh$Country=="Bosnia-Herzegovina", "black", "slategray3"), lwd=1)
grid(NA, NA, col="#333333", lty="dotted", lwd=0.2)
axis(1, tick=FALSE, cex.axis=0.9)
axis(2, tick=TRUE, cex.axis=0.9)
GDP <- df.fh$GDPChange1989
mod1 <- lm(df.fh$TypeAVol~GDP)
abline(mod1, col="gray2", lwd=3)
newx <- seq(min(df.fh$GDPChange1989), max(df.fh$GDPChange1989), 0.01)
CI <- predict(mod1, newdata = data.frame(GDP=newx), interval = c("confidence"), level = 0.95)
x <- sort(newx)
line1 <- sort(CI[,2], decreasing=TRUE)
line2 <- sort(CI[,3], decreasing=TRUE)
lines(x, line1, type="l", lty=2, lwd=.5, col="red")
lines(x, line2, type="l", lty=2, lwd=.5, col="red")
text(df.fh$GDPChange1989[df.fh$Country=="Bosnia-Herzegovina"], df.fh$TypeAVol[df.fh$Country =="Bosnia-Herzegovina"], df.fh$pclabel[df.fh$Country =="Bosnia-Herzegovina"], cex=0.8, pos=2)